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We have studied the effect of nonzero bulk viscosity with peak near the lattice QCD predicted 
crossover temperature Tco ~ 175MeV on charged particle transverse momentum spectra and elliptic 
flow. The Israel-Stewart theory of 2nd order causal dissipative relativistic fluid dynamics is used 
to simulate the space time evolution of the matter formed in Au-Au collisions at ^/sNN=200 GeV 
assuming longitudinal boost invariance. A systematic comparison of temperature, transverse veloc- 
ity, spatial and momentum anisotropy evolution of the ideal, bulk and shear viscous fluid has been 
carried out. Two different temperature dependent forms of ("/s and a constant 77/s was used. Both 
the bulk and shear viscous correction to the ideal freezeout distribution function are included. The 
dissipative correction to the freezeout distribution for bulk viscosity was calculated using Grad's 
fourteen moment method. From our simulation we show that the method is applicable only for 
C/s < 0.005 X tj/sIkss for freezeout temperatures 130 and 160 MeV. 

PACS numbers: 12.38.Mh ,47.75.-|-f, 25.75.Ld 



I. INTRODUCTION 

One of the most interesting results coming from RHIC 
heavy ion program is the observation that hot QCD mat- 
ter created in Au-Au collisions behaves like an almost 
ideal fluid [1-4]. 

Relativistic hydrodynamics has been a successful the- 
ory to describe the bulk properties of the QCD medium 
produced in nuclear collision. Hydrodynamic simulation 
of nuclear collision at RHIC indicates that the shear vis- 
cosity of QCD plasma is very small. However the ex- 
tracted value of shear viscosity largely depend on the ini- 
tial condition [6]. The theoretical estimation of various 
kinetic coefficient for a QCD plasma becomes enormously 
complex due to strong coupling. String theoretical calcu- 
lation based on the ADS/CFT correspondence predicts 
a lower limit on the ratio of shear viscosity to entropy 
density as 77/5 > l/in [7]. 

The bulk viscosity ^ can in general be of the same order 
of magnitude as shear viscosity rj. But until recently, the 
bulk viscosity was neglected in the study of fluid dynam- 
ics of QCD matter. There are important special cases 
in which C is very much smaller, or vanishes altogether 
[8] . In [8] author shows that a simple gas of structureless 
point particles will have negligible bulk viscosity in the 
extreme-relativistic or non relativistic limits. However, it 
must be stressed that a vanishing bulk viscosity is the ex- 
ception, rather than the rule, for general imperfect fluids 
[8],[9]. 

Existence of non zero bulk viscosity in QCD plasma 
near the critical temperature was reported in [10-12] 
from lattice QCD calculation. The theoretical estima- 
tion for transport coefficients in QCD plasma can also 
be found in [13-18]. In [13] a lower and upper bound 
of shear viscosity was given and bulk viscosity associ- 
ated with the plasma to hadron transition was estimated 
in the relaxation- time approximation. Several sources 
of shear and bulk viscosity was discussed in [18] with 
the emphasis that the bulk viscosity is associated with 
the chemical nonequilibrium. The effect of bulk viscos- 



ity on freezeout and HBT puzzle was studied in [19]. 
In spite of all these theoretical calculations, there are 
several model dependency in the estimate of the trans- 
port coefficients. One can introduce the dissipative ef- 
fects in hydrodynamics by treating transport coefficient 
like shear and bulk viscosity as input parameter and ob- 
tain their values in phenomenological study by compar- 
ing to the experimental data. In the following, we will 
consider effect of bulk viscosity on hydrodynamic evolu- 
tion and subsequent particle production. The temper- 
ature dependence of C/s is different for leading order 
pQCD and AdS/CFT calculation [10]. In pQCD cal- 
culation C = 1577(T) (1/3- c2(r))^ and in AdS/CFT 
C ^ ??(r) (1/3 — c^(r)), where Cs is the speed of sound 
in the medium, 77 is the shear viscosity and T is the tem- 
perature. C/s(T) is different among the available lattice 
calculations [10-12]. In this study we will describe in 
detail the results obtained from the recently extended 
version of our code by considering two different tempera- 
ture dependent forms of C/s and a constant rj/s. In addi- 
tion to shear viscous correction to freezeout distribution 
function we have also considered the bulk viscous cor- 
rection to the ideal freezeout distribution function. The 
bulk viscous correction dfbuik was calculated by using 
Grad's 14-moment method for a multicomponent gas as 
discussed in [20] . Considering a small value of bulk stress 
H at the freezeout it was shown in [20] that bulk viscosity 
has a non-negligible effect on particle spectra and ellip- 
tic flow coefficient. Thus it is important to consider the 
bulk viscous evolution and the dissipative correction to 
the ideal freezeout distribution function. We will con- 
centrate in this work to study the effect of bulk viscosity 
on pt spectra and elliptic flow of charged hadron with 
for two different freezeout temperature T/=130 and 160 
MeV. We want to remind that the calculation of Sfbuik 
for a multicomponent hadron gas is non trivial, for a re- 
cent calculation of the dissipative correction 5f to the 
single particle distribution function in leading-log QCD 
and in several simplified model see [21]. 
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The paper is organized as follows, in section II, we dis- 
cuss the formalism of bulk viscous hydrodynamics in the 
context of Israel-Stewart 2nd order theory. Section III 
deals with the input (jy/s, C/s(T), initial condition, re- 
laxation time, equation of state and freezeout condition) 
required for the 2-l-lD viscous hydrodynamics simulation. 
In section IV we present the results of our simulation 
which includes the time evolution of bulk viscous stress 
Il{x,y) as well as its spatial average {Il{x,y)), and the 
effect of C/s on evolution and observables. A comparative 
study between ideal, shear and bulk viscous fluid evolu- 
tion has also been done. The effect of the bulk viscous 
correction to the ideal freeze-out distribution function on 
charged pion's elliptic flow and pr spectra is discussed. 
Finally in section V we present a summary of our study. 

II. FORMULATION 

The most widely studied theory of relativistic causal 
dissipative hydrodynamics is due to Israel-Stewart (I- 
S) [22-24]. In the present work we follow the I-S for- 
malism of viscous hydrodynamics in 2-l-lD to study the 
effect of bulk viscosity on experimental observables. For 
completeness we start with the brief description of the 
formalism which is followed by a detail discussion on im- 
plementing bulk viscosity in the fluid evolution as well as 
the corresponding dissipative correction to the freezeout 
distribution function. 

For a simple fluid in equilibrium or in adiabatic limit, 
the energy momentum tensor T'^'' and particle current 
Ni^ has the following forms [8] 



particle current and energy-momentum tensor can be 
written as, 

AN^"' = -n-^ (4) 

e + p ^ ' 

AT^"" = U{u''u'' - g"") + . (5) 

In Eq.4, is the heat conduction current, 11 is the 
bulk viscous stress and w^'^ is the shear stress tensor. The 

space time evolution of the fluid is governed by the energy 
momentum and particle number conservation laws, 

a^T"^ = 0, (6) 
d^Ni" = 0. (7) 

Apart from these two conservation laws, the second 
law of thermodynamics restrict the entropy four vector 
in such a way that the rate of entropy production 
per unit volume should always be positive for all possible 
fluid configurations, 

5^5'' > (8) 

For an imperfect fluid the entropy four vector can 
be decomposed into two parts. 

S-^ = S^g + AS" (9) 

where S^^ is the equilibrium part and AS" is the dissi- 
pative correction. The equilibrium entropy four current 
is given by, 



eg — 

N" = 



(e + p)u"u'' 
nu" 



(1) 
(2) 



where e ,p and n are total energy density, pressure, 
and particle number density, u" is the velocity four- 
vector, normalized so that w/'w,^ = 1. Choice of hydro- 
dynamic velocity is arbitrary. There are two commonly 
used comoving frame to fix the u, (i) Landau-Lifshitz 
frame where u" is parallel to the energy four flow and 
(ii)Eckart frame where the u" is parallel to the A^^. We 
are using the Landau-Lifshitz frame, which is more ap- 
propriate than the Eckart frame, in baryon free central 
rapidity region. 

Due to the presence of dissipative processes the energy 
momentum tensor and the particle four flow of ideal fluid 
gets modified as follows: 



rp^v rp^v _|_ ^rp^v 

m = m„ + Am 



(3) 



AT'"' and AN^^ are dissipative corrections to the en- 
ergy momentum tensor and particle 4-current. In the 
Landau-Lifshitz frame, the dissipative correction to the 



S" = 



pu'" 



(10) 



If we neglect the heat conduction and consider bulk 

and shear viscosity as the only dissipation mechanism, 
the non-equilibrium correction to the entropy four cur- 
rent can be written as, 



AS^ = -^u"U^-^u"n^,n' 



(11) 



where the co-efficients and /32 are function of energy 
density and number density. Their exact values can be 
determined from kinetic theory [22]. Later we will iden- 
tify that these co-efficients are related to the relaxation 
time of bulk and shear viscosity respectively. 

The AS''' is constructed in terms of various orders of 
the gradient of u". For a first order theory the second 
law of thermodynamics d^S" > can be satisfied by 
postulating the following form of bulk 11 and shear w"" 
stress 



n 



ail' 



-CO 



<''u"> 



(12) 
(13) 



3 



where C and rj are the positive transport coefScients, bulk 
and shear viscosity respectively. 9 is the expansion rate 
to be defined later. 

It can be shown that the first order theory violates 
causality [25]. For example if, in a given fluid cell, 
thermodynamic forces vanish, corresponding dissipative 
fluxes also vanishes instantly all over the volume. Causal- 
ity violation of dissipative hydrodynamics is corrected in 
2nd order theories [22]. According to this theory the 
shear and bulk viscous stress follows the relaxation equa- 
tions. 



0.10 



(14) 
(15) 



where D — u^d^ is the convective time derivative, 9 = 
d.u is the expansion scalar and V^^u^^ = ^(V'^m'^ -|- 
y^u^) — ^{d.u){g^'^ — u^^u^) is a symmetric traceless ten- 
sor. One can identify C/3o and 277/32 respectively with the 
relaxation time for bulk and shear stress tensor, Tn — CPo 
and = 277/32. 

The relaxation equation 14,15 for bulk and shear vis- 
cosity have additional terms according to the kinetic the- 
ory calculation [26] and has the following form. 



dh 



1 



[2f]V 



TT 



1 _ „ 1 

rn^ ■ 2 



(16) 



[n + cv^7i^ + -CTna^(^)].(i7) 



Assuming longitudinal boost invariance, earlier, we 
have solved the relaxation equations for shear stress ten- 
sors, using the code 'AZHYDRO-KOLKATA'. Detafls of 
the code can be found in [24] . We have extended the code 
to include the bulk viscosity. Relaxation equations for 
bulk and shear viscosity are solved simultaneously with 
the energy-momentum conservation equations. Details 
of the equations solved are given in the appendix A. 



III. EQUATION OF STATE, VISCOSITY 
COEFFICIENTS AND INITIAL CONDITIONS 
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FIG. 1: (Color online)Two different form of temperature de- 
pendence of (/s. (a)Form-l: (^/s in the QGP phase (T 
>175 MeV) is calculated by using PQCD formula C/s = 
15f (r)(l/3 — Cs(r))^, where was calculated from recent 
lattice data [32]. In the hadronic phase (T < 175 MeV)C/s 
is parametrized from [35]. (b)Form-2: This form is taken 
from [27], where in the QGP phase ("/s was obtained from a 
different lattice calculation [12] and C/s in hadronic phase is 
from [35]. Red dashed line is the KSS bound [7] of the shear 
viscosity to entropy density ratio rj/s ~ l/47r. 



calculation [32] . Entropy density of the two phases were 
smoothly joined at T = Tc = 175MeV by a step like 
function. 

The thermodynamic variable pressure, energy density 
etc. are then calculated by using the standard thermo- 
dynamic relations. 



A. Equation of State 

Equation of state is one of the important input to hy- 
drodynamic model. Through this input, the macroscopic 
hydrodynamic models make contact with the microscopic 
world. In the present simulations we have used an equa- 
tion of state with cross-over transition at Tc — 175MeV, 
developed earlier [42] . The low temperature phase of the 
EoS is modeled by hadronic resonance gas, containing all 
the resonances with Mres <2.5 GeV. The high tempera- 
ture phase is a parametrization of the recent lattice QCD 



p{T) = r s{T')dT' 

JQ 

e{T) = TS{T)~P{T) 
B. Bulk and shear viscosity coefficients 



(18) 
(19) 



As discussed in the introduction the exact form of 
the ('/•s(T) is quite uncertain. In this work we choose 
two different temperature dependent form. The form-1 
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FIG. 3: (Color online)The spatially averaged bulk viscous 
stress n(a;, y) in the transverse plane as a function of evo- 
lution time for two different initialization of n(a;,y)and (a) 
for ^/s(T)form-l (b)<^/s(T)form-2. The simulation with zero 
initial bulk stress Tl{x, y) = is marked as blue. Sim- 
ulation for NS initialization Il{x,y) = —(,8 is marked as 
red. The solid,dashed dotted and dashed lines represent the 
simulation with different relaxation time for n(a;, y) namely 
rn=1.0,0.5and 0.1 times r^. 



FIG. 2: (Color online)Top Panel: the transverse profile of 
bulk stress 'n{x,y) = —C^6 at initial time to=0.6 fm/c. Bot- 
tom Panel: Il{x,y) at a later time r — 11.7 fm. 



is shown in Fig. 1(a). It is constructed in the following 
way, for the QGP phase we use the formula derived in 
pQCD calculation C/s = 15f (T) (l/3 - c2(r))^ The 
squared sound speed c^iT) is calculated using the rela- 
tion cl{T) = dp{T)/de{T) |^ where the pressure p(r) and 
energy density e(T) are obtained form the lattice QCD 
calculation [32]. The rj/s{T) used in this calculation is 
same as given in [33] which was obtained using lattice 
data [34]. The C/s for the hadronic phase used in form-1 
is the same as calculated in [35]. Where it was calcu- 
lated by using hadron resonance gas model including all 
known hadrons and their resonances up to mass 2 GeV 
with finite volume correction to thermodynamic quanti- 
ties. It also included an exponentially increasing density 
of Hagedorn states in the mass range of 2-80 GeV. 
The form-2 is shown in Fig. 1(b). This form was taken 



from the reference [27]. For the QGP phase the C/s was 
taken from a different lattice calculation [12]. For the 
hadronic phase the Q/ s was taken from [35]. 

The peak value of form-2 is ^10 times larger than the 
peak value in form-1. Though both form of Q/ s shows 
a peak near the crossover temperature {T^o ^175 MeV), 
their dependency on temperature is slightly different in 
the QGP phase. The red dashed line in both Fig. 1(a) and 
(b) shows the KSS bound of shear viscosity to entropy 
density ratio [7]. 

For comparison purpose, we will also show some sim- 
ulation result with only shear viscosity. However, in this 
demonstrative simulation, we have neglected the tem- 
perature dependence of shear viscosity to entropy ratio 
and perform the simulations with the AdS / CFT minimal 
value, r]/s = I/A-k. 
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C. Relaxation Time 

Solution of relaxation equations for bulk and shear 

stress tensors requires to specify the relaxation time (rn) 
for the bulk stress and (t^-) for the shear stress tensor. In 
principle, relaxation times rn and could be calculated 
from the underlying kinetic theory, which for strongly 
coupled QCD plasma, is a complex problem. Relax- 
ation times Tn and Ttt was calculated in [22] for sim- 
ple relativistic Boltzmann, Bose and Fermi gases with 
mass m using Grad 14 moment approximation in rela- 
tivistic kinetic theory. For a Boltzmann gas, in the non- 
relativistic limit (/3 = ^ oo), m = (Po ~ zT^p 
and = 2?7/32 « v/P^ In the photon limit (/3 — >■ 0), 
rn = C/3o - C^^i^^y and = 277^2 = ^. Note that 
in the photon limit, the mass term appear in the the de- 
nominator with a quadratic power. The relaxation time 
become very large. For very large relaxation time, bulk 
stress will evolve very slowly with time. To circumvent 
the problem, for bulk stress also, one generally use the 
relaxation time for the shear stress tensor. In our simu- 
lation the Tn is either a constant or same as for the shear 
viscosity. 



D. Initial Conditions 



where A^pa,.((b, x, y) and Ncou{h,x,y) are the transverse 
profile of participant numbers and binary collision num- 
bers respectively, cq in Eq.20 docs not depend on the 
impact parameter of the collision. It corresponds to the 
central energy density in b=0 impact parameter colli- 
sion. Generally eo is fixed to reproduce the experimental 
charged hadron multiplicity or the pt spectra in a central 
collision. Analysis of experimental data in -^5=200 GeV 
Au+Au collisions indicate that the energy density of the 
central fluid is cq «30 GeV/ fm^. In the present study, 
we fix eo = SOGeV/ fm^. x in Eq.20 is the hard scat- 
tering fraction. It was shown in [36] that with x=0.13 
the experimental ccntrality dependence of charged mul- 
tiplicity data can be best explained at RHIC energy. The 
effect of varying 'x' has been studied in [38], [39]. For the 
present study, we fix x = 0.13. For the present study 
we carry out simulation for Au-Au collision at an impact 
parameter b=7.4 fm unless stated otherwise. 

In this study wc will use Navier-Stokes(NS) initializa- 
tion of bulk stress n(a;, y) = —Q9 unless stated otherwise. 
One can also initialize n(a;, y) by assuming a zero value 
at the initial time t^. Figure 2 shows the initial Y\{x.y) 
in the transverse plane with the Navier-Stokes initializa- 
tion. 



The initial conditions used here includes the ini- 
tial energy density profile (e(.T,y)) in the transverse 
plane, the initial time (r^) ,the transverse velocity profile 
{vx{x,y),Vy{x,y)), shear and bulk stresses in the trans- 
verse plane {'K{x,y),Yi{x,y)) respectively at Tj. The val- 
ues of these parameter are given in the table below. One 
also need to specify the freezeout conditions to stop the 
hydrodynamics evolution, this will be discussed in a later 
section. 



TABLE I: Initial conditions for 2-l-lD viscous hydrodynamics 
calculation. 



Parameters 


Values 


£0 


30 (GeV//m^) 


n 


0.6 fm 


Vx{x,y),Vy{x,y) 


0.0 




2??/3Ti, 


n{x,y) 


-CO 



The initial energy density profile in the transverse 
plane was parameterized in a two component Glauber 
model. At an impact parameter b, the transverse energy 
density is obtained from the following two component 
form 



e(b,x,y) = eo[(l - x)NpaTt{^,x,y) + xNcon{h,x,y)] 

(20) 



E. Preezeout 



In the following, we assume that the fluid undergoes ki- 
netic freeze-out at a fixed temperature Tp- We have con- 
sidered two freeze-out temperature, Tf=130 MeV and 
160 MeV. In hydrodynamical model, one generally as- 
sumes that prior to kinetic freeze-out, fluid undergoes 
a 'chemical freeze-out' at temperature Tchem > Tp, be- 
yond which the particle ratios remain unaltered. In the 
present simulations, we have used an lattice QCD based 
EoS. Lattice QCD calculations arc for zero baryon den- 
sity fluid. We therefore assume that the chemical equi- 
libration is maintained throughout the evolution. The 
proton-antiproton ratio will not be correctly reproduced 
in the model. However, protons constitute only 5% of 
the total charged particle yield. For other particles e.g. 
pion, kaon, the model simulation can correctly reproduce 
experimental numbers. Moreover, in the present simula- 
tions, we have studied the effect of bulk viscosity on hy- 
drodynamic evolution and particle production. We have 
not attempted fits to experimental data. 

We have used the Cooper-Frye algorithm [43] to eval- 
uate particle spectra from the freeze-out surface. We 
have included the dissipative correction to particle spec- 
tra. For detail see the appendix B. 
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FIG. 4: (Color onlme)The rate of cooling of the fluid for ideal, 
shear and bulk viscous evolution is shown here for three dif- 
ferent location ,r=0 fm( solid line), at r=2.8 fm (dashed line) 
and at r=5.6 fm (dashed dotted line). Red line is for ideal 
fluid evolution, blue and pink lines are for bulk and shear 
viscous evolution respectively. Top panel (a) the calculation 
with (/s form-1. Bottom panel(b) with ("/s form-2. 



IV. RESULTS 
A. Effect of bulk viscosity on fluid evolution 

The space-time evolution of Tl{x, y) is governed by the 
relaxation Eq. 17. The relaxation time in Eq. 17 controls 
how fast the stress n(x, y) relaxes to its instantaneous 
equilibrium value. The top panel of figure 2 shows the ini- 
tial n(a;, y) with Navier-Stokes initialization in b=7.4 fm 
Au-|-Au collisions for form-1 of bulk viscosity. In b=7.4 
fm collision, collision zone is asymmetric. The bottom 
panel of the same figure shows the n(a;, y) at a later time 
T=11.76 fm. The relaxation time is rn = Ttt. Within a 
span of '^lO fm, bulk stress is decreased by a factor of 
100. One also note that with time, anisotropicity of the 
bulk stress is also decreased. Pressure gradient is more 
along the minor axis (cc-direction in Fig. 2), and fiuid ex- 
pands more in that direction reducing the anisotropicity. 
Form-II of C/s also give similar results. 

In Fig. 3(a) we have shown the temporal evolution of 
spatially averaged bulk stress, (n(a;, y)) for two differ- 




FIG. 5: (Color online) Time evolution of spatially averaged 
transverse velocity {vt) for ideal and viscous simulation. The 
red solid line is for ideal fluid evolution and the dashed dot and 
dotted lines are for bulk and shear viscous evolution respec- 
tively. Top panel (a) is the simulation with form-1. Bottom 
panel (b) calculation with C,/s form-2. 



ent initializations (i) Ti{x, y)—Q at the initial time Ti(blue 
lines) and (ii) Navier-Stokes(NS) initialization n(a;, y) = 
—C,0 (red lines) at t^. Results are shown for three differ- 
ent relaxation time, rn = 1.0 (solid line) ,0.5 (dashed dot 
line) and 0.1 (dashed lines) times r^. Tj^=3ri/2p is the re- 
laxation time for shear viscous stress estimated for a rel- 
ativistic Boltzmann gas[22]. Late time evolution of bulk 
stress hardly depend on the initialization. Even when the 
bulk stress is initialized to zero value, it quickly reaches 
the Navier-Stokes value. Time by which the bulk stress 
reach the Navier-Stokes value depend on the relaxation 
time. The shorter relaxation time tjj means the system 
reaches its equilibrium Navier-Stokes value faster, and 
in the limit of rn — Eq.l7 transforms to relativistic 
Navier-Stokes equation (Eq.l2). From Fig. 3 one can see 
that (n(a;, y)) with zero initialization takes least time to 
attain its instantaneous equilibrium value (red line) for 
the smallest value of rn. In Fig. 3(b), same results are 
shown for the form-2 of C/s. Results are similar to that 
obtained for form-1 of C/s. 

The rate of cooling of the fiuid element for two different 
form of C/s at three different location in the reaction zone 
for ideal (red line), bulk (blue line) and shear (pink line) 
viscous evolution are shown in Fig 4(a) and (b). Fig 4(a) 
is the simulation for C/s form-1 and 4(b) is for (^/s form- 
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FIG. 6: (Color online)The evolution of spatial anisotropy with 
time. Solid red line is the ideal evolution, dashed dotted line 
is for bulk viscous and dashed line indicates the shear viscous 
evolution, (a) Simulation with ("/s form-1. (b)Simulation 
with ("/s form-2 
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FIG. 7: (Color online)The evolution of momentum anisotropy 
with time. The lines shows here are the same as explained in 
figure 6. Fig.(a) is the result for form-1 of C,l s and (b) is for 
form-2 of C^j s. 



2. The rate of cooling is different at various points in the 
reaction zone. There is no noticeable change in the rate 
of cooling due to bulk viscosity compared to ideal fluid 
evolution. For evolution with ry/s = 0.08 decreases the 
rate of cooling in the early time for the central region. 
The difference between the effect of bulk and shear stress 
on rate of cooling is due to the difference in magnitude 
of n and -K^^ . In the peripheral region the temperature 
variation is almost same for the ideal, shear and bulk 
viscous evolution. Early in the evolution (2-3 fm) the 
fluid expansion is mainly in the longitudinal direction 
and follows Bjorken cooling law T^t — constant. At a 
later time the transverse expansion leads to a different 
slope for the cooling rate. 

Fig 5(a) and (b) shows the temporal evolution of 
the spatially averaged transverse velocity (((wt))) of 
the fluid with form-1 and form-2 of C/s respectively. 
The space averaged transverse velocity is defined as 

{{vt)) = '^'^^^("^1)"" ■ Sohd red line is for ideal fluid 
and the dashed dotted and dotted line is for bulk and 
shear viscous evolution respectively. Because of the 
reduced pressure gradient in the bulk viscous evolution 
compared to ideal fluid, the corresponding ((wt)) gets 
reduced. The reduction in ((wt)) at later time (after ~' 
8 fm) is more for the C/s form-2 compared to form-1. 



The shear viscosity on the other hand increases the 
pressure gradient in the transverse direction and reduced 
the longitudinal pressure at the early time of evolution. 
Because of the enhanced pressure gradient ,the ((wt)) 
for shear viscous evolution is increased compared to 
ideal fluid. 

In a non-central collision between two identical nuclei, 
the collision zone is non-spherical. The spatial eccentric- 
ity Ex is defined as 



ex = 



(21) 



is a measure of spatial deformation of the fireball from 
spherical shape. A zero value of means the system 
is spherical, < Ea; < 1 indicates elliptic shape with 
major axis along Y direction, and e^; < means the ma- 
jor axis along X direction. The angular bracket ((...)) 
implies an energy density weighted average. For b=7.4 
fm collision, the evolution of Ex with time (t) is depicted 
in Fig. 6(a) and (b). The solid red line corresponds to 
the temporal evolution of for ideal fluid, the dashed 
dotted and dashed line are for bulk and shear viscous 
fluid evolution respectively. Because of the reduced pres- 
sure gradient in bulk viscous evolution, the initial spatial 
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FIG. 8: (Color online)Negative charged pion's pr spectra for 
Au+Au collision at impact parameter h—YA fm.The solid red 
line is the spectra from ideal fluid. The dashed dot and dotted 
lines are for fluid with bulk viscosity with form-f and form- 
2. The inset figure shows the ratio of correction to the pr 
spectra due to bulk viscosity to ideal evolution. 



FIG. 9: (Color online) Negative charged pion's V2 for Au+Au 
collision at impact parameter b=7.4 fm.The solid red line is 
the «2 from ideal fluid. The dashed dot and dotted lines are 
for fluid with bulk viscosity with form-f and form-2. The 
inset figure shows the ratio of correction to the V2 due to bulk 
viscosity to ideal evolution. 



deformation {sx ~ 0.28) takes a longer time to change 
its shape for the bulk viscous evolution compared to the 
ideal fluid evolution. At a later time the change in the Sx 
is more pronounced for form-2 of C/s compared to form- 
1. Shear viscosity does the opposite to the transverse 
expansion, initially it enhances the transverse velocity 
and the spatial deformation Ex reduces at a much higher 
rate compared to the ideal fluid evolution. 

Similar to the spatial anisotropy, one can define the 
asymmetry of fireball in momentum space. The momen- 
tum space anisotropy Sp is defined as 



The simulated elliptic flow V2 in hydrodynamic model 
is directly related to the temporal evolution of the 
momentum anisotropy. In fact in ideal hydrodynamics 
V2 oc £p. The evolution of Sp as a function of time 
is shown in Fig. 7(a) and (b) for C,/s form-1 and Q/s 
form-2 respectively. Compared to ideal evolution (solid 
red line), the bulk viscous evolution (dashed dot line) 
results in a reduced value of momentum anisotropy 
around freezeout time(^ 12 fm/c). The change in tp for 
bulk viscous evolution compared to ideal fluid occurs 
after t ~ 3-4 fm/c. Around this time most regions of 
the fluid element reaches the temperature range ^175 
MeV where C/s has maximum value. 



B. Effect of bulk viscosity on particle spectra and 
elliptic flow 

1. Without freeze- out correction 

We first discuss the change in particle spectra and el- 
liptic flow due to bulk viscosity. Non-equilibrium cor- 
rection to equilibrium distribution function is neglected. 
The simulated charged pion pT spectra is shown in the 
Fig. 8. The red solid line in 8 is the result obtained in 
ideal fluid simulation, the dashed dot and dotted line is 
obtained with form-1 and form-2 of bulk viscosity respec- 
tively. Due to the smaller transverse flow in presence of 
bulk viscosity we get a steeper pT spectra compared to 
the ideal fluid evolution. However, the change is small. 
The inset figure shows the relative correction in the pt 
spectra SN/Neq due to the bulk viscosity compared to 
ideal simulation. Where SN= N^uik — N^q- The cor- 
rection to the Pt spectra is small, at px ~2 GeV the 
correction is ^ 10% with form-2 of Q/s. It is even less 
for form-1. The simulated elliptic flow V2 of tt" pro- 
duced in Au-Au collision at impact parameter b=7.4 fm 
collision are shown in Fig. 9. The red solid line is the 
result for ideal hydrodynamics. The dashed dotted and 
dotted lines are V2 with bulk viscosity with form-1 and 
form-2 respectively. Here, again, we have not included 
the non-equilibrium correction to the distribution func- 
tion. The inset of 9 shows the relative change in the 
V2 with bulk viscosity compared to ideal fluid evolution. 
The relative correction in V2 due to the form-1 of C,/ s is 
within ^ 4% for the pT range 0-3 GeV and for form-2 the 
relative correction is within ~ 8%. It appear that if the 
non-equilibrium correction to the equilibrium distribu- 
tion function is neglected, both the form of bulk viscos- 
ity introduces relatively small correction to the particle 
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P, (GeV) 



FIG. 10: (Color online) The charged pion pr spectra for im- 
pact parameter b=7.4fm Au-Au collision. Ideal fluid evolu- 
tion (black solid line) and various values of C^/s form-1. (a) 
Freezeout temperature T/^ISO MeV (b) Tf = 160 MeV. 



spectra and elliptic flow. 



2. With freeze-out correction 

It has already been discussed that there are two kinds 
of dissipative correction to the ideal fluid simulation. 
First the energy momentum tensor contains a viscous 
correction and the freezeout distribution function is also 
modified in the presence of dissipative processes. So far 
in this paper all the bulk viscous simulation results are 
obtained for dissipative correction to the energy momen- 
tum tensor only. In this section we discuss the correction 
to the freezeout distribution function. We have employed 
Grad's fourteen-moment methods for freezeout dissipa- 
tive correction as described in [20]. The implementation 
of this method to our 2+VD viscous code '"AZHYDRO- 
KOLKATA"' is briefly discussed in the appendix B. 

Figure 10 (a) and (b) shows the pT spectra of pion for 
freezeout temperature T/=130 and 160 MeV for four dif- 
ferent values of C,/ s form-1. The black solid line in Fig. 10 
is the pt spectra for ideal fluid evolution in 20-30% Au- 
Au collision, the other lines are for bulk viscous evolu- 
tion with varying values of form-1 C,/s. The red dashed 
line is the simulated spectra with form-1 C/s, whereas 
long dashed, dashed dotted and dotted lines are results 
for 0.1,0.05 and 0.01 times the form-1 C/s. It appears 



_\ \ 


(a) T,=130MeV 
C/s 

0.1 ^/S 

0.05 ^/s 

0.01 C/s 








\\ , 7 





(b) T,=160MeV 

C/s 

o.U/s 

0.05 11% 

0.01 lj% 


'///////// 









0.0 0.5 1.0 1.5 2.0 2.5 3.0 

Pt (GeV) 

FIG. 11: (Color online)The relative correction to the pt spec- 
tra due to the bulk viscosity compared to ideal fluid evolu- 
tion. Dissipative correction in both evolution as well in freeze- 
out distribution function has been included. Shaded region 
shows the relative correction of 50%. (a) Freezeout tempera- 
ture r/=130 MeV (b) r/=160 MeV. 



that for QCD motivated bulk viscosity over entropy ra- 
tio, (f)QCD = 157(^ — Cg)^, Grad'sl4-moment method 
introduces very large correction to the spectra. At low 
Pt, compared to ideal fluid, vr" yield is increased by a 
factor of 10 or more. Corrections are comparatively less 
if bulk viscosity is reduced. For very small bulk viscos- 
ity, 0.01 times form-1 C/s, the spectra is very close to the 
ideal fluid evolution in the pt range < < ^GeV. 

To put things in perspective for the present calcu- 
lations relative to the existing calculations in [20], we 
should compare the relative change in invariant yield of 
negative pions versus px between ideal simulation and 
bulk viscous simulation for both cases. In order to do 
that, we have to compare with the corresponding calcu- 
lations done in the present work with input bulk viscosity 
to entropy density of 0.01 x C/s at Tj = 160 MeV. Such 
a value of C/s is chosen so as to have a similar magnitude 
of n over the freeze-out surface as used in [20] . We find 
that the relative corrections are similar. 

The specific form of the dissipative correction to the 
ideal freezeout distribution function considered here (see 
appendix B) leads to a large negative correction to the 
Pt spectra for large values of px- Depending on the value 
of C/s the dissipative correction due to the bulk viscos- 
ity results in a negative invariant yield above a certain 
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value oi pt- A negative value of particle number is un- 
physical, we will omit the px range in the subsequent 
plots from where the particle number becomes negative. 
As discussed earlier, freeze-out correction is obtained un- 
der the assumption that the non-equilibrium correction 
to the distribution function is small than the equilib- 
rium distribution function, 5 f hulk « feq- It is then 
implied that the relative correction (SN/Neq) is small for 
Israel-Stewart's hydrodynamics to be applicable. Figure 
11 (a) and (b) shows the corresponding relative correc- 
tion {SN/Neq) to the pt spectra due to the bulk viscosity. 
The shaded band in the figure corresponds to the rela- 
tive correction of 50%. We consider here a correction of 
magnitude greater than 50% indicates the breakdown of 
the the freezeout correction procedure. From Fig. 11 one 
can see that the pt spectra changes drastically in bulk 
viscous evolution (solid black line) with C/s form-1. Only 
for viscous simulation with bulk viscosity to entropy den- 
sity ratio less than 0.01 times (/s the relative correction 
is less than 50%. The results impose a severe constraint 
on the bulk viscosity, | < 0.01(|) qcd. 

We find from Fig. 10 and 11 that the value of 5 f bulk 
is greater for the simulations with Tf = 130 MeV than 
those for Tf = 160 MeV. This can be understood in 
the following manner. As shown in the last equation 
of Appendix B both the values of 11 on the freeze-out 
surface as well as the values of the prefactors Dq, Bq 
and -Bq at a given value of Tf determines the magni- 
tude of 5 f bulk- The average magnitude of 11 decreases 
from about -5 x IQ-^GeV/fm^ at Tf = 160 MeV to 
-2 X lO-^GeV/fm? at Tf = 130 MeV. However the pref- 
actor values as given in Table II increases with decrease in 
temperature. So the observed Tf dependence of Sfbuik is 
due to the interplay of both the temperature dependence 
of n and the prefactors. 

Fig. 12 shows the elliptic flow of pion for 20-30% col- 
lision centrality as a function of pr for ideal and bulk 
viscous evolution. The lines represents the same condi- 
tions as described in Fig. 10. The relative correction to 
the V2 is defined in the same way as for the pr spec- 
tra is shown in Fig. 13 (a) and (b) for T/=130MeV and 
160 MeV respectively. The V2 puts a more stricter con- 
straint than Pt spectra on the allowed value of C/s for 
the applicability of Grad's 14 moment method. We find 
the relative correction within 50% for bulk viscosity to 
entropy density ratio less than 0.05 times C/s form-1. 

In all the above calculation, for the relaxation time of 
bulk viscosity, we have used the Boltzmann estimate for 
the relaxation time for shear viscosity rn = TV . We have 
investigated the effect of relaxation time on the pT spec- 
tra and elliptic flow. For relaxation time rn=0.1, 1 and 
5 times Ttt, we have solved the hydrodynamic evolution 
and computed 7r~ pr spectra and elliptic flow. Results 
are shown in Fig. 14. If relaxation time is decreased by a 
factor for 10 from rn = r^r to rn = 0.1 Ttt, spectra or el- 
liptic flow hardly changed. If relaxation time is increased 
by a factor of 5, at large pr yield is reduced marginally. 
Effect is more pronounced on elliptic flow. V2 is in- 




(b)T,=160MeV I I j 




Pt (GeV) 



FIG. 12: (Color online) Same as figure 10 but for elliptic fiow 

creased. Increased fiow with increasing relaxation time 
can be understood as follows: bulk stress evolve compar- 
atively slowly with increased Tn and the non-equilibrium 
correction at the freeze-out is increased. Increased non- 
equilibrium correction will lead to increased V2- 

V. SUMMARY 

The effect of bulk viscosity on pion px spectra and el- 
liptic fiow was studied by numerically solving 2-)- Id rela- 
tivistic viscous hydrodynamics equations. Two different 
parametrize form of C/s(T) was used along with constant 
Tj/s. To construct C/s in the partonic phase we use lat- 
tice data, (/s in hadronic phase is calculated using a 
hadron resonance gas model including hagedorn states 
with limiting mass of 80 GeV. 

A comparative study of ideal, shear and bulk viscous 
evolution is also done. We observe that the time varia- 
tion of the temperature of the fiuid remains similar for 
ideal and bulk viscous evolution. Whereas in presence of 
shear viscosity the cooling rate is reduced. Because of the 
reduced pressure due to finite C/s ,the transverse veloc- 
ity slightly decreases compared to the ideal fluid around 
freezeout time. The shear viscosity on the other hand in- 
creases the transverse pressure which results in a higher 
transverse velocity compared to the ideal evolution. The 
observable consequence of the above fact is reflected in 
the slope of the px spectra of pion. The time evolution of 
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FIG. 13: (Color online) Same as figure 11 but for elliptic flow 

V2- 



spatial eccentricity is almost unchanged in ideal and bulk 
viscous evolution. Due to the large transverse velocity 
in shear viscous evolution the spatial deformation shows 
a rapid change compared to ideal fluid. The momentum 
anisotropy in the shear viscous evolution grows at a much 
slower rate compared to ideal evolution. The change in 
the time variation of momentum anisotropy in the bulk 
viscous fluid with respect to ideal case is observed at a 
late time of the evolution. This will have consequences 
to the observed elliptic flow of the hadrons. 

The modification of the freezeout distribution function 
according to Grad's moment method has been consid- 
ered. The change in pT spectra and V2 in bulk viscous 
evolution with respect to ideal simulation with no cor- 
rection to the freezeout distribution is within 5-10% de- 
pending on the form of C/s(T). However a large correc- 
tion to both the px spectra and elliptic flow was observed 
for bulk viscous simulation with dissipative correction to 
the freezeout distribution function. Combined study of 
Pt spectra and V2 of pion in 20-30% Au-Au collision with 
full bulk viscous evolution puts a constraints on the ap- 
plicability of Grad's 14 moment method. We find the 
relative correction within 50% for bulk viscosity to en- 
tropy density ratio less than 0.01 times C/s form-1. 
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Appendix A: Evolution equations 

The energy momentum tensor including shear and bulk 
viscosity has the following component in 2+1 dimension 
in r, a;, j/, r\ co-ordinate. 



rj-17 



(e + p + n)7i-(p + n) 

T^^ = {e+p + U)jlv, 

= {e + p + n)-flvy- 



TT 



(Al) 



The conservation equations df^T^" — in 2+1 dimension 
are written in the following forms to solve them numer- 
ically. The three energy-momentum conservation equa- 
tion along with four relaxation equation (1 eq. for bulk 
and 3 eq. for shear viscosity) for viscous stress tensor are 
solved by using the '"SHASTA" '[46] algorithm. 



H + rV") 



(A2) 



drT^- + d^iT^^v,) + dyir^'vy) = 

-d,{p + n + n--- v^n--) - dyiW-' - n^'-'vy) (A3) 
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where 



Where T^"' =TT^"' a.ndp = TP ,U = tU,v^ = ^, = _a 1a™ 6» 



3 
1 



^77 J'TT • 2 



3 



The relaxation equation for bulk and shear tensor takes a^^ = —^i^xU^^ dyU^ — u^Du^ — u^Du^] 

the following form ^ 

— A^^6» (A6) 

an an an 3 ^ ' 



■J-[n + c^ + inrnax + CTnDf^)]. (as) 



and 



dxa: I a i Q ^ { ^xx '^^„xx\ ^ rxx 

a^TT™ + v^axTT™ + %aj,7r™ = ^- (tt™ - 2r?c7™) - -7i™ -fr= 2u2' [7r°^£)Mo - 7r^^£>«.; - 7r™£'%] (A7) 

T,r7 7 

^rTT^^ + i^cS^tt"*' + z;j,aj,7r"^ = ^ (7r"« - 277a"«) - i/f^ 

7V7 7 

I 



Here D = is the convective time derivative and 

Q = dp,u^ is the expansion scalar, r^r is the relaxation 
time for shear stress, = 2ry/32 and rn = C/^o is the 
relaxation time for bulk viscous stress. 

We solve only above three component of shear stress. 
The dependent shear stress tensor components can be 
obtained from the independent ones. Using the prop- 
erties that (i) T^^^ is transverse to and (ii) tt''^ is 
traceless, g^^i^^v = 0,the dependent shear stress tensor 
components can be obtained as, 

TT^^ = ?;^7r^^ + VyTT'^y (A9) 

■K-^y = v^-K^'y + vy-K^y (Alo) 
TT^^ = vI-k'"'' + vl^yy + iv^vy-n'^y (au) 

+2v^Vyn''y (A12) 

Appendix B: Bulk viscous correction to freezeout 
distribution 

In Cooper Frey prescription, the particle momentum 
distribution is obtained by integrating the single particle 
distribution function over the freezeout hyper surface S^. 



where E is the energy, g is degeneracy, and is four 
momentum. In ideal hydrodynamics the fluid is in local 
thermal equilibrium and the distribution function is the 
one particle equilibrium distribution function feq{x,p), 

f{x,p) = r{x,p) = ^3 ^ 1 (B2) 

with inverse temperature /3 = l/T and chemical poten- 
tial /i. g is the degeneracy factor. The (±) arc respec- 
tively for fermions and bosons. For dissipative fluids, the 
system is not in local thermal equilibrium. In a highly 
non-equilibrium system, distribution function f{x,p) is 
unknown. If the system is slightly off-equilibrium, then 
it is possible to calculate correction to equilibrium dis- 
tribution function due to (small) non-equilibrium effects. 
Slightly off-equilibrium distribution function can be ap- 
proximated as, 

f{x,p)=feg{x,p)+Sf (B3) 

Where Sf = Sfbuik + Sf shear « f represents the dissi- 
pative correction to the equilibrium distribution function 
/eg, due to bulk viscosity and shear viscosity. 

There are different methods available to calculate 
the dissipative correction to the distribution function 
[20, 21, 27, 37, 41, 44]. In order that the energy- 
momentum tensor remains continuous across the freeze- 
out surface the functional form of the df must be such 
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that the Landau matching condition should be satisfied; 
u^l^T^'^Ui, — [41]. For bulk viscosity, the dissipative 
correction for a multicomponent system was calculated 
by using Grad's fourteen moment method in [20]. Fol- 
lowing [20] the dissipative correction for bulk viscosity 
Sfbulk can be written in the following way, 



Sfbulk = -feq{^ + £feq) X 



n 



Where the prefactor Do , Bo and Bq are temperature 
dependent constant. Here we have dropped the index 
'i' of particle species for simplicity. The Landau match- 
ing condition is satisfied with the present form of correc- 
tion to the ideal distribution function. In [20] the prcfac- 
tors Dq, Bo and Bq was calculated for a multicomponent 
hadron gas. We use their estimated values (given in the 
table II) in this calculation for two different freezeout 
temperature T/ = 130Mey and T/ = 160MeV. 

TABLE II: Prefactors for two different temperature 





DoiGeV-^) 


Bo{GeV-^) 




130 MeV 


9.10x10* 


1.12x10' 


-3.27x10* 


160 MeV 


2.01x10* 


1.66x10* 


-7.84X 10^ 



Shear viscous correction to equilibrium distribution 



function is well known, 



shear — feqi^ ~t~ ^feq) 



2{e + p)T^ 



PllPuT^' 



(B4) 



The correction to the ideal spectra "^^i^^y^ due to the 
bulk viscosity is 



= / di:,p^5f,.,,[p^u,. T) (B5) 

Now the four momentum of the fluid element is p^ = 

{mTCOshy,Px,Py,mTsinhy) where ttit = \J vn^ + p\ and 
the momentum rapidity is y = ^In^z^- The freeze-out 
hypersurface in (r, x, y, jf) co-ordinate is 

rfSju = ^^TCOshr],—^^,—^^mTsinhri^ Tfdxdydr) 

and 

p^.dTif^ = (rnTCOsh{r] — y)— Pt-Vtt/ j Tfdxdydrj 

Using these relationship into equation B5 and after 
some algebra we have the final form of the correction 
to the invariant yield due to the bulk viscosity which is 
given as 



J 



d'^Pxdy \ A 4 2 2 

-Pt-Vt-t/ |^fco(n/3_L) + ^fc2(n/3j_) + 62^1 (n^_L) + hkQ{nl3s_) 



where j3±_ — niTjP 



A = 2^(Tl)"+^e 



(^Bq — Bo^ 



nj3{-iVTPT) 



rriT [Oaj + {p^'v^ +p"vy) Bo - 2Boj^ {p^'v^ + p^Vy)^ 
Boml - Doi (p'^v, + p^vy) - BopI (l + j^l) - BopI (l + i^l) 
- 2Bop''vxpyvyj^ + Boj^ {p^v^f + 2Boj^p^Vxpyvy + {p'^Vyf 



14 



K. Adcox et al. [PHENIX Collaboration], Nucl. Phys. A [25 
757, 184 (2005). 
J. Adams et al. [STAR Collaboration], Nucl. Phys. A [26 
757, 102 (2005). 
I. Arscnc et al. [BRAHMS Collaboration], Nucl. Phys. A [27 
757, 1 (2005). 

M. Gyulassy and L. McLerran, Nucl. Phys. A 750, 30 [28 
(2005). [29 
A. K. Chaudhuri, Phys. Lett. B681, 418-422 (2009). 
T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and [30 
Y. Nara, Phys. Lett. B 636, 299 (2006). [31 
P. Kovtun, D. T. Son and A. O. Starinets, Phys. Rev. 
Lett. 94, 111601 (2005). [32 
S. Weinberg, Astrophys. J. 168, 175 (1971). [33 
Landau L.D,and Lifshitz E.M. 1959. Fluid Mechanics, 
page 308-312. [34 
H. B. Meyer, Phys. Rev. Lett. 100, 162001 (2008). 
S. Sakai and A. Nakamura, PoSLAT 2007, 221 (2007). [35 

F. Karsch, D. Kharzeev and K. Tuchin, Phys. Lett. B 
663, 217 (2008). [36 
P. Danielewicz, M. Gyulassy, Phys. Rev. D31, 53-62 [37 
(1985). [38| 
J. I. Kapusta, arXiv:0809.3746 [nucl-th]. 

S. Jeon, Phys. Rev. D 52, 3591 (1995). [39 
A. Wiranata and M. Prakash, Nucl. Phys. A 830, 219C 
(2009). [40 
M. Prakash, M. PraJiash, R. Venugopalan and G. Welke, 
Phys. Rept. 227, 321 (1993). [41 
K. Paech and S. Pratt, Phys. Rev. C 74, 014901 (2006). 

G. Torrieri, B. Tomasik and I. Mishustin, Phys. Rev. C [42 
77, 034903 (2008). 
A. Monnai and T. Hirano, Phys. Rev. C 80, 054906 [43 
(2009). 

K. Dusling and T. Schaefer, arXiv:1109.5181 [hep-ph]. [44 
W. Israel, Nonstationary Irreversible Thermodynamics: 
A Causal Rolativistic Theory, Annals of Physics 100,310- [45 
331 (1976). [46^ 
A. Muronga, Phys. Rev. C 76, 014909 (2007). 
A. K. Chaudhuri, arXiv:0801.3180 [nucl-th]. 



W. A. Hiscock and L. Lindblom, Phys. Rev. D 31, 725 
(1985). 

E. Molnar, H. Niemi and D. H. Rischke, Eur. Phys. J. C 
65, 615 (2010). 

G. S. Denicol, T. Kodama, T. Koide and Ph. Mota, Phys. 
Rev. C 80, 064901 (2009). 

D. Kharzeev and K. Tuchin, ,JHEP 0809, 093 (2008). 
J. W. Chen and J. Wang, Phys. Rev. C 79, 044913 

(2009) . 

D. Davesne, Phys. Rev. C53, 3069-3084 (1996). 

H. Song, "Causal Viscous Hydrodynamics for Relativistic 
Heavy Ion Collisions," [arXiv:0908.3656 [nucl-th]]. 

S. Borsanyi et al, JEEP 1011, 077 (2010). 

H. Niemi, G. S. Denicol, P. Huovineu, E. Molnar, 

D. H. Rischke, Phys. Rev. Lett. 106, 212302 (2011). 

A. Nakamura, S. Sakai, Phys. Rev. Lett. 94, 072305 

(2005). 

J. Noronha-Hostler, J. Noronha and C. Greiner, Phys. 
Rev. Lett. 103, 172302 (2009). 

T. Hirano and Y. Nara, Phys. Rev. C 79, 064904 (2009). 

P. Bozek, arXiv: 11 10.6742 [nucl-th]. 

V. Roy and A. K. Chaudhuri, Phys. Rev. C 81, 067901 

(2010) . 

P. F. Kolb, U. W. Heinz, P. Huovinen, K. J. Eskola and 
K. Tuominen, Nucl. Phys. A 696, 197 (2001). 
X. G. Huang, T. Kodama, T. Koide and D. H. Rischke, 
Phys. Rev. G 83, 024906 (2011). 

M. Luzum and J. -Y. Ollitrault, Phys. Rev. C 82, 014906 

(2010) . 

V. Roy and A. K. Chaudhuri, Phys. Lett. B 703, 313 

(2011) . 

F. Cooper, G. Fryc and E. Schonborg, Phys. Rev. Dll, 
192 (1975). 

K. Dusling and D. Teaney, Phys. Rev. C 77, 034905 
(2008). 

A. K. Chaudhuri, J. Phys. G G37, 075011 (2010). 

J. P. Boris and D. L. Book, J. Comput. Phys. 11, 38 

(1973). 



